4  Using Meta-regressions to Compare CDD across Species or Sites

4.1 Overview

Following from Section 4. How does CDD vary across species, abiotic gradients or in time? in the main text, here, we demonstrate one approach to comparing the strength of CDD across species using a meta-analysis framework. The same approach can be used to compare CDD across sites, plots, or any other unit of interest, as long as it is possible to generate reliable estimates of the strength of CDD (suitable sample sizes, etc.). As noted in the main text, “Correct propagation of uncertainty in CDD estimates requires meta-regressions in frequentist or Bayesian frameworks. Weighted regressions (e.g., lm, lmer, gam) can also estimate how CDD varies e.g., across latitude or between species with different life-history strategies, but incorrectly estimate the associated uncertainty.”

We use a subset of the BCI seedling data (Comita et al. 2023) of only 30 species to compare how species abundance is related to the strength of CDD at the species level. This is the same dataset used in section 2. This analysis is for demonstration purposes only, and biological conclusions should not be made about the results, given this is only a small subset of the data.

An advantage of using meta-regressions over simple weighted regressions is that the models are able to simultaneously account for uncertainty in species-specific CDD estimates as well as systematic differences in species’ CDD when regressed against a predictor via the inclusion of random effects. A simple weighted regression, on the other hand, assumes that there is only one error for which the relative strength is known when regressing the estimates against a predictor. The latter can lead to incorrect weighting of species in the metaregression.

In this tutorial, we use relative average marginal effect (𝑟𝐴𝑀𝐸) calculated in the previous chapter as our response variable, calculated separately for each species. 𝑟𝐴𝑀𝐸 in this case estimates the relative increase in the probability of annual mortality with the addition of one new conspecific neighbor, while keeping total densities constant. Positive numbers indicate a relative increase in mortality with an increase in conspecific density, a signature of NCDD. In principle, any metric of the strength of CDD can be used, though care must be taken to ensure that the metrics are comparable across species and sites (see main text for more information).

We use the popular metafor package to fit the meta-regression models.

Note

The following code is adapted from the latitudinalCNDD repository by Lisa Hülsmann.

4.2 Load libraries and data

Code

# Load libraries
library(tidyr)
library(dplyr)
library(readr)
library(ggplot2)
library(here)
library(metafor)

# Load in species abundances for BCI data subset
abundances <- read_csv(
here("./data/BCI seedling data - 30 species - abundance 2023_05_18.csv")
  )

# Load marginal effects calculations from previous section
load(here("./data/mortality.Rdata"))

# Subset down to just equilibrium change
rAME <- rAME %>%
  filter(change == "equilibrium")

# Join marginal effects and abundance data
rAME <- left_join(rAME, abundances, by = c("sp" = "spp"))

# Add in average abundance for 'rare' species
rAME$abundance[rAME$sp == "data_deficient_seedling"] <- 
  mean(abundances$abundance[abundances$spp %in% 
                              nsp$spp[nsp$data_deficient]])

# Log transform abundance to use in models
rAME$log_abundance <- log(rAME$abundance)

---

Let’s take a quick look at the data set we’ll be working with:

Code

head(rAME, n = 10)

       term    estimate   std.error offset change.value      change     sp
1  con_dens  0.12190548 0.122096554      1          + 1 equilibrium ACALDI
2  con_dens  0.46332386 0.140120748      1          + 1 equilibrium AEGICE
3  con_dens  0.02026319 0.003624488      1          + 1 equilibrium BEILPE
4  con_dens  0.42931069 0.087277501      1          + 1 equilibrium CAPPFR
5  con_dens  0.04121409 0.042571758      1          + 1 equilibrium CECRIN
6  con_dens  0.10440629 0.241547627      1          + 1 equilibrium CORDLA
7  con_dens  0.11658356 0.200294935      1          + 1 equilibrium DAVINI
8  con_dens  0.01113154 0.123539647      1          + 1 equilibrium DESMAX
9  con_dens  0.23375216 0.227252222      1          + 1 equilibrium HEISCO
10 con_dens -0.06498043 0.061086207      1          + 1 equilibrium JUSTGR
   abundance log_abundance
1        169      5.129899
2        153      5.030438
3       5574      8.625868
4       1604      7.380256
5         51      3.931826
6        248      5.513429
7         27      3.295837
8         55      4.007333
9         56      4.025352
10       136      4.912655

---

4.3 Reformat data for model fitting

First, we use the ‘escalc’ function in the metafor package to essentially repackage our data frame into a format used in the meta-regression model fitting. Since we already calculated our effect size (𝑟𝐴𝑀𝐸), we just pass through the 𝑟𝐴𝑀𝐸 estimate and corresponding standard error using the ‘GEN’ option for the ‘measure’ argument, rather than calculating an effect size within the ‘escalc’ function.

Code

# Reformat data for model fitting
# Set measure to generic, which passes the observed effect sizes or 
# outcomes via the yi argument and the corresponding sampling 
# variances via the vi argument (or the standard errors via the sei 
# argument) to the function.
    dat_meta = metafor::escalc(measure = "GEN", 
                               yi = estimate, # observed outcomes
                               sei = std.error, # standard errors
                               slab = sp, # label for species
                               data = rAME)

---

4.4 Fit meta-regression model

Next, we use the ‘rma’ function to fit a meta-regression model, where 𝑟𝐴𝑀𝐸 is our response variable (renamed as yi in the previous step) and log species abundance is our predictor. While not shown here, it is possible to fit mixed effects meta-regression models with the ‘rma.mv’ function. We suggest consulting the extensive documentation for the metafor package for further details.

Code

# Fit model
metamod = metafor::rma(yi = yi,
                       vi = vi,
                       mods = ~ log_abundance,
                       method = "REML",
                       data = dat_meta)

4.5 Print model summary

Code

summary(metamod)

Mixed-Effects Model (k = 22; tau^2 estimator: REML)

  logLik  deviance       AIC       BIC      AICc   
  2.2364   -4.4727    1.5273    4.5145    3.0273   

tau^2 (estimated amount of residual heterogeneity):     0.0166 (SE = 0.0090)
tau (square root of estimated tau^2 value):             0.1287
I^2 (residual heterogeneity / unaccounted variability): 91.39%
H^2 (unaccounted variability / sampling variability):   11.62
R^2 (amount of heterogeneity accounted for):            0.00%

Test for Residual Heterogeneity:
QE(df = 20) = 74.0420, p-val < .0001

Test of Moderators (coefficient 2):
QM(df = 1) = 0.3751, p-val = 0.5402

Model Results:

               estimate      se     zval    pval    ci.lb   ci.ub    
intrcpt          0.1935  0.1401   1.3815  0.1671  -0.0810  0.4681    
log_abundance   -0.0141  0.0231  -0.6125  0.5402  -0.0594  0.0311    

---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

4.6 Plot the estimated of rAME for all species with a forest plot

In this case, a forest plot shows the estimates of the strength of CDD for individual species, here ordered by least to most abundant going from top to bottom.

Code

forest(metamod, 
       header = "Species", 
       xlab = "rAME",
       order = log_abundance)

4.7 Model diagnostics

The plot method displays model diagnostics of the meta-regression model (more info here) in addition to a forest plot.

Code

plot(metamod)

4.8 Model predictions of how species abundance is related to strength of CDD

Here, we generate predictions and corresponding confidence intervals for how our predictor of interest is related to the related to strength of CDD using the ‘predict’ function. The y-axis here indicates the relative increase in annual mortality probablity with the addition of one conspecific neighbor. Higher values indicate stronger negative conspecific density dependence. In this example, our predictor of interest is species abundance. We also scale the size of the points based on their weight in the meta-analysis, with larger points indicating higher weights. Note that this analysis is for demonstration purposes only, and biological conclusions should not be made about the results, given this is only a small subset of the data.

Code

# Generate a prediction dataframe
pred <- expand_grid(log_abundance = seq(min(dat_meta$log_abundance, 
                                            na.rm = TRUE), 
                                        max(dat_meta$log_abundance, 
                                            na.rm = TRUE),
                                        length.out = 50))

pred$abundance <- exp(pred$log_abundance) # Back transform abundance

# Bind predictions to dataframe
pred <- cbind(pred, predict(object = metamod, 
                            newmods = pred$log_abundance))

# Extract observed values
observed_values <- broom::augment(metamod)

# Add in variance estimates to be able to scale size of points by 
# amount of variance in estimate
observed_values <- left_join(observed_values,
                             dat_meta %>% 
                               dplyr::select(sp, vi, log_abundance),
                             by = c(".rownames" = "sp"))



# Set abundance values for x axis
abundances_x_axis <- c(25, 50, 100, 1000, 5000)

# Plot prediction
ggplot(pred, aes(x = log_abundance, y = pred)) + 
  geom_ribbon(aes(ymin = ci.lb, ymax = ci.ub), 
              fill = "steelblue2", alpha = 0.75) + 
  geom_line() + 
  geom_hline(yintercept = 0, lty = 2) + 
  labs(x = "Species Abundance", y = "rAME of CDD") + 
  scale_x_continuous(breaks = log(abundances_x_axis),
                     labels = abundances_x_axis) + 
  # Add observed points
  geom_point(data = observed_values,
             aes(x = log_abundance, y = .observed, size = 1/vi), 
             alpha = 0.75) + 
  theme_bw(15) +
  theme(legend.position = "none")

In principle, the approach outlined in this appendix of estimating marginal effects followed by a meta-regression analysis can be used to test for a correlation between strength of CDD and other species-level traits (e.g., wood density, leaf area, etc) or to determine how site-level CDD varies with site-level variables, such as latitude (Hülsmann et al. 202x), precipitation, soil fertility, etc. 

Comita, Liza S, Salomón Aguilar, Stephen P Hubbell, and Rolando Pérez. 2023. “Long-Term Seedling and Small Sapling Census Data from the b Arro c Olorado i Sland 50 Ha f Orest d Ynamics p Lot, p Anama.” Ecology 104 (9): e4140.